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Synchronization processes in populations of identical networked oscillators are in the focus of intense studies 
in physical, biological, technological and social systems. Here we analyze the stability of the synchronization of 
a network of oscillators coupled through different variables. Under the assumption of an equal topology of con¬ 
nections for all variables, the master stability function formalism allows assessing and quantifying the stability 
properties of the synchronization manifold when the coupling is transferred from one variable to another. We 
report on the existence of an optimal coupling transference that maximizes the stability of the synchronous state 
in a network of Rossler-like oscillators. Finally, we design an experimental implementation (using nonlinear 
electronic circuits) which grounds the robustness of the theoretical predictions against parameter mismatches, 
as well as against intrinsic noise of the system. 

PACS numbers: 89.75.Fb,89.75.Hc,05.45.-a 


I. INTRODUCTION 

Synchronization processes on complex networks has re¬ 
ceived a lot of attention during the last decades IBEl. The 
interplay between the dynamical evolution of oscillators and 
their local interactions (as given by the complex topology of 
a network) usually results in non-trivial phenomena of rele¬ 
vance to physical, biological, technological and social sys¬ 
tems. First introduced by Pecora and Carroll (H, the Master 
Stability Function (MSF) is nowadays one of the main theoret¬ 
ical methods for the study of network synchronization. MSF 
is indeed a powerful tool to analyze the stability of the syn¬ 
chronization manifold when identical systems of oscillators 
are diffusively coupled. Originally applied to undirected net¬ 
works, the MSF approach has been later extended to inves¬ 
tigate enhancements and optimization of complete synchro¬ 
nization in weighted and asymmetric topologies (see dEl, 
and references therein). 

In 0 the authors stated the so-called heterogeneity para¬ 
dox, i.e. the fact that heterogeneous networks, wherein dis¬ 
tances between nodes are relatively short, are less stable, in 
terms of synchronization, than their homogeneous counter¬ 
parts. Soon after, a proper and adequate weighting of the link 
strengths was shown to overcome this paradox, based again on 
concepts sparkling from the MSF formalism 017]. Follow¬ 
ing works, have shown how different network’s topological 
features influence the stability of the synchronous state, such 
as: heterogeneity of the node degree, degree-degree correla¬ 
tions, average shortest-path, betweenness centrality or cluster¬ 
ing. These latter studies indicate that altering the structure of 
a network may result in maximizing the stability of the syn¬ 
chronous state, thus achieving a maximally stable synchro¬ 


nization structure 10. Enhancement of the networks’ syn- 
chronizability can also be achieved by the application of ge¬ 
netic algorithms increasing the stability of the synchronized 
state. In this case, the networks self-organize by disconnect¬ 
ing the hubs and connecting peripheral nodes, thus increasing 
the homogeneity and leading to what is known as entangled 
networks El. 

In our study, we report the enhancement of the stability of 
complete synchronization of an ensemble of dynamical units, 
when coupled simultaneously in different dimensions. We are 
concerned with a multivariable coupling, where the dynamical 
systems are coupled through different dimensions according 
to a certain structure of connections (see Fig for a schematic 
illustration). In particular, we consider a generic dynami¬ 
cal system whose associated vector state x (with x G 
evolves according to x = f (x). Each one of the m state vari¬ 
ables of the dynamical system at a given node can be coupled 
to the corresponding variable of any of the other systems (i.e., 
nodes) of the network. 

Equivalently, we can think of our system as a network with 
I < m layers, each one accounting for the structure of cou¬ 
plings at each variable of the system. This multilayer point 
of view illustrated in Figj^is, in fact, just accounting for a 
multivariable coupling between the nodes of a network, nev¬ 
ertheless it will help us to provide a more concrete represen¬ 
tation of the structure of the system, and possible connections 
to applications. So we will make use of it at certain points. If 
the coupling between oscillators does not include some of the 
state variables, i.e. I < m, the topology of the corresponding 
layers to those variables would be trivially given by a zero ad¬ 
jacency matrix, so we would not consider them to be proper 
layers (as is the case of the layer corresponding to variable z 
in Fig[^. For simplicity, we consider a bidirectional coupling 
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FIG. 1: (Color online) Network of iV = 6 dynamical systems 
coupled through different variables. In order to better observe the 
coupling introduced at each variable, each node i is split into differ¬ 
ent layers, each one corresponding to a variable of the system (X, 
Y and Z), whose complete dynamical state is given by the vector 
X = {x,y,z) obtained from the combination of the states of its vari¬ 
ables at each layer. Note that, in this particular example, only vari¬ 
ables X and y are used to couple the systems and, in turn, that the 
topology of the coupling network is the same for both variables. 


between the same variables of each system (i.e. each layer is 
an undirected network). This is illustrated in Fig. [^with an 
example of the case I = 2 and m = 3. 

Interestingly, our framework connects with the so-called 
hypernetwork formalism introduced by Sorrentino fiOl . In 
this latter work, the author shows that a MSF approach to hy¬ 
pernetwork synchronization is possible when the Laplacian 
matrices m of different layers (accounting for the coupling 
through each variable) have the same basis of eigenvectors, 
i.e., when they are simultaneously diagonalizable. This is a 
condition that has been shown to be fulfilled for two layers 
in three cases: (i) the Laplacian matrices of the different lay¬ 
ers are commuting (a condition that automatically allows for a 
MSF approach whatever the number of layers if the Laplacian 
matrices form a pairwise commuting set), (ii) one of the two 
layers is unweighted and fully connected, or (iii) one of the 
two layers has an adjacency matrix of the form Aij = aj with 
i,j = 1,...,N. Additionally, Ref. ||T2| contains an extension 
of the approach in cni to more general topologies by mak¬ 
ing use of a simultaneous block diagonalization of Laplacian 
matrices corresponding to different layers, thereby decreasing 
the dimensionality of the linear stability problem. 

In our work, we consider the topology to be the same in 
each layer, trivially falling, from the of view of hypernet¬ 
works, into category (i) of Ref. cni. This way, we present 
a study on how the stability of the synchronous state is en¬ 
hanced by finding an optimal balance for the coupling be¬ 
tween the different variables in a network of identical oscilla¬ 
tors with multivariable coupling. On the one hand, we provide 
results based on extensive numerical simulations of networks 
of Rossler-like oscillators to show the applicability of the pro¬ 
posed ideas and how the MSF can help us to find the adequate 
balance between the couplings that optimizes the stability of 
the synchronous state of a network. On the other hand, by 


constructing an electronic version of the model, we show that 
these predictions are in good agreement with the experimental 
evidences in spite of the idealizations used in the theoretical 
treatment. 


II. COUPLING THROUGH DIFFERENT VARIABLES 
WITH IDENTICAL TOPOLOGY 


In this section we explain how stability of the synchronous 
state can be enhanced by engineering a multivariable coupling 
function between nodes in a network and what balance be¬ 
tween coupling variables is the most adequate. For the sake 
of concreteness we focus on a set of Rossler-like oscillators 
ca coupled to their neighbors through both the x and the y 
variables, whose dynamics evolve according to the following 
equations: 


N 

Xi = -ai{xiY Pyi^Vzi- {I - dr)(Jx^'^afj[xj - Xi]) 

i=i 

N 

Vi = -Ol 2 {-lXi + {l-5)yi- dr(JY(t>YYj [Vj - Vi]) 

i=i 

{^Xi H“ 5 (1) 

where ai = 500, 0^2 = 200, 0^3 = 10,000, /3 = 10, T = 20, 

7 = 50, ^ = 8.772, /i = 15, = 20, 0 = 50 and (1 - dr)(Jx 

and drdY account for the coupling strengths of variables x 
and y. As we explain below, this chaotic oscillator has the 
highly non-trivial characteristic of being quite robust when 
implemented in electronic circuits. The adjacency matrices 
and contain the topology of each of two layers, each 
one accounting for the coupling through the x and x variables. 
Elements afj and are one when nodes i and j are con¬ 
nected and zero otherwise. With these parameters the oscil¬ 
lators display chaotic dynamics due to the nonlinearity intro¬ 
duced in Ga,., which consists on a piecewise function defined 
as: 




r 0 Xi <3 

\ p{xi- 3 ) Xi> 3 


( 2 ) 


The coupling between oscillators is here controlled by two 
parameters: a being the coupling strength and dr controlling 
how the coupling strength is distributed between variables x 
and y. This way, dr = ^ {dr = Y) leads to a coupling re¬ 
stricted to variable x (y), while a sweep of dr in the interval 
[ 0 , 1 ] allows for a weighted combination of both x and y vari¬ 
ables. Notice that the role of the parameter dr is therefore that 
of exploring the consequences of unevenly distributed cou¬ 
pling on the stability of synchronization, which in the past 
were the object of specific studies in space extended systems 
d and weighted monolayer graphs ca Here however, 
the difference in the weight assigned to each variable intro¬ 
duced by dr implies a different balance of two system’s vari¬ 
ables in the coupling configuration, but does not affect the 
un-directionality of each one of the network’s links. 

We now apply the MSF formalism to study how the amount 
of coupling effectively distributed among the two coupling 
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variables affects the stability of the synchronous state of the hx,y : ^ , the dynamics of each node is then given 

network. Denoting the coupling functions of each variable as by: 


Xi 


N N 

f(Xi) + {l-dr)crx'^Afj[hx{Xj) -hx{Xi)]+ dray'y'^Alj [hy(xj) - hy(Xj)] 
i=i i=i 


JV / 

N 

\ N / 

N 

f(xi) + (1 - dr)ax E] 


hx(Xj) + dray J2\AF- Sij 


.=1 V 

_i=i 

j i=i V 

_i=i 


N N 

f(Xi) - {1 - dr)crx'^C^jhx{Xj) - dray'^CjjhYiXj) 
i=i i=i 


hr(xj) 


(3) 


X Y 

where Sij stands for the Kronecker delta, and C-' are the 
Laplacian matrices m describing the coupling through vari¬ 
ables X and y respectively. If we consider that ay = crx = cr 
and we restrict our analysis to the case of identical coupling 
topologies for all variables of the system, i.e. = A?^ = 
A, and, in turn, = £, Eq. [^reads: 

N 

Xi=f(xi) - a^£ijh(xj) (4) 

i=i 

where the coupling function is h(xj) = (1 — d^)hx(xj) — 
drily (xj). This way, Eq. |^is basically the classical equa¬ 
tion describing the evolution of a diffusively coupled systems, 
with the particularity that the coupling function h depends on 
the parameter dr accounting for how the total amount of cou¬ 
pling is divided between the coupling variables of the system. 
The dependence on dr leads to a parametric MSE that de¬ 
scribes the stability of the synchronization manifold. Varying 
the value of dr we obtain a family of MSEs that allows to eval¬ 
uate how the stability of the synchronized manifold is affected 
by shifting the coupling from one layer to the other. 

The independent variable u of the MSE is related with the 
eigenvalues of the Laplacian matrix {a = and the syn¬ 
chronization manifold will be stable when, for all eigenval¬ 
ues of the Laplacian matrix 7 /., the corresponding a leads to 
a value of the MSE that is negative (for a given value of cr). 
Taking into account that 71 = 0 (since the Laplacian is a zero 
row sum matrix) and that 72 < 73 < ... < 7Ar, dynamical 
systems can be classified, depending on the shape of its corre¬ 
sponding MSE. The classification includes: i) class I systems, 
whose MSE is always positive for any value of u so that the 
system cannot be synchronized for any topology or coupling 
strength, ii) class II systems, when the MSE is positive for low 
values of u and becomes negative when a threshold value uj 
is achieved (being the synchronization manifold stable when 
(7^2 > III’ when the MSE has two zeros at 

TAi and z^//, leading to a stable synchronized manifold when 
ui < <772 and uii > cf^n simultaneously. 

Now we investigate the synchronization properties of the 


system depending on the distribution of the coupling strengths 
among the coupling variables. In Eig. [^a) we show the MSE 
obtained for the Rossler-like system coupled through the x 
variable (dr = 0 , black solid line), the y variable (dr = 1 . 0 , 
black dashed line) and simultaneously the x and y variables 
(colored lines). A range of dr values is swept in order to show 
the gradual changes in the stability of the synchronous state as 
the coupling of one or the other variable is enhanced. As ex¬ 
pected, coupling introduced only through variable x leads the 
system to be class III while it becomes class II when coupled 
only through the y variable CSIIIII. Interestingly, when the 
coupling is introduced simultaneously through the two vari¬ 
ables we obtain a MSE that is not a linear combination of the 
isolated layers. In particular, the sweep of dr leads to a family 
of MSEs that are class II (at least for the ten cases plotted in 
Eig.J^a)), thus synchronizing the network when the condition 
crj 2 ^ is fulfilled, which is in principle achievable for any 
topology by just using a sufficiently high a. 


Figure [^a) also points out an important consequence, 
namely that using a multivariable coupling (while maintain¬ 
ing the overall coupling strength) leads to values of the MSE 
that are more negative than those obtained when coupling the 
systems through a unique variable, indicating that the stabil¬ 
ity of the synchronized manifold is higher for multivariable 
couplings. This is somehow to be expected, since a simple 
change of coordinates would lead the dynamics of the Rossler 
system to be described as a combination of the actual x, y and 
^ variables. However, by studying how the MSE modifies its 
shape as a function of the combination of the variables of the 
system we can find the regions with the highest stability. This 
fact can be observed in Eig. [^b), where the MSE is plotted 
as a function of dr for three different values of u . We can 
observe a minimum of the MSE at intermediate values of dr 
for a = 2 and a = 3, i.e., sl dr value where the stability of 
the synchronized manifold is the largest. This result is also 
observed within a range of values of u (not shown here). 
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FIG. 2: (Color online) Master Stability Function MSF(z/) as a func¬ 
tion of the coupling fraction dr between the variables to be coupled. 
In (a), solid and dashed black lines are, respectively, the MSF for a 
coupling through variables x and y exclusively. Colored lines show 
the MSF obtained for different values of dr, as indicated in the leg¬ 
end. In (b), MSFs for three different values of the parameter 
z/ = 1, z/ = 2 and z/ = 3. Note that, for a sizable range of dr, 
the multivariable coupling leads to lower values of the MSF, indi¬ 
cating a region where the stability against external perturbations is 
higher. 


III. NUMERICAL SIMULATIONS 

With the aim of testing the predictions of the MSF, we sim¬ 
ulate a network of A/" = 6 bidirectionally coupled Rossler-like 
systems, see Fig. and Eq. Q- We introduce an additive 
Gaussian noise term (with zero mean) in the x variable whose 
strength is given by the parameter r]. Next, we calculate the 
synchronization error e of the network as a function of the 
coupling a and the noise strength r], for different values of dr. 
This way, we are able to check the robustness of the stability 
of the synchronized manifold depending on how the coupling 
is distributed among the x and y variables of the system. To 
facilitate the comparison with the MSF predictions, we plot 
the results as a function of u = ( 772 , instead of a, being 


(a) 
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FIG. 3: (Color online) Robustness of the synchronized state. Syn¬ 
chronization error e as a function of the noise strength 77 and the 
coupling (rescaled as z/ = ^ 772 ). Panels correspond to: (a) dr = 0, 
(b) dr = 1 and (c) dr = 0.5. The network coupled through both the 
X and y variables (dr = 0.5) is the one showing better performance 
against noise perturbations. The synchronization error in the system 
is computed as '^i<j \xi~XjV where the normalizing fac¬ 

tor corresponds to the total number of oscillator pairs in the network. 
Left plots correspond to the MSF obtained for each value of dr. 


72 = 2.382 for the particular network structure shown in Fig. 
[2 This choice can be justified on the grounds that all the im¬ 
portant information from the point of view of the synchroniza¬ 
tion is given by the MSF value of the eigenmode associated to 
72 , as the eigenmode associated to 76 is never pushed beyond 
the boundaries of the synchronization region within the ex¬ 
perimentally accessible range of coupling strengths. Figure 
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FIG. 4: (Color online) Experimental setup. On the left we show a schematic representation of the coupling topology of the 6-circuit network. 
The coupling is adjusted using one digital potentiometers X9C104, whose parameters Cu/d (Up/Down resistance) and Cstep (increment of the 
resistance at each step) are controlled by a digital signal coming from a DAQ Card, PO.O-PO. respectively. The outputs of the circuit are sent 
to a set of voltage followers that act as a buffer and, then, sent to the analog ports ( AI 0 ; AI 1; ... ; AI 5) of the same DAQ Card. The ports 
DAC generate the 6 noisy signals to test the robustness of the network. The whole experiment is controlled from a PC with Labview Software. 


shows the two-dimensional plots e(z/, 77 ) obtained for dr = 0 
(coupling only through the x variable), dr = 0.5 (x and y cou¬ 
plings are equally active) and dr = 1 (coupling only through 
the y variable). For each of these three cases, the coupling 
strength a is modified and the synchronization error between 
all oscillators is calculated (see caption of Figj^for details). 
We consider that the network is out of synchrony for synchro¬ 
nization errors e > 0.6 (red regions in Fig. |^. 

When the coupling is introduced through the x variable 
(Fig. [^(a)), the Rossler-like oscillators behave as a class III 
for any value of the noise strength (within the range 0 < 77 < 
5), i.e. the system only synchronizes for intermediate val¬ 
ues of the coupling strength u. When the noise strength 77 
is increased, the synchronization error increases, leading to 
a complete loss of the synchronization for large values of 77 . 
For couplings such that u 1.2 or larger, the network be¬ 
comes unstable in the sense that the strong coupling makes 
the individual oscillators abandon the basin of the attractor, 
and their dynamics blow up. This type of instability is to be 
expected, as a similar phenomenon is observed in individual 
Rossler-like oscillators for some initial conditions (see, e.g., 
Ref. msi). The reader should notice that no such problems 
beset the computation of the MSF, as the maximum Lyapunov 
exponent transverse to the synchronization manifold is com¬ 
puted from effectively one individual oscillator along a tra¬ 


jectory that follows the attractor dynamics. Thus, there is no 
contradiction in the fact that the MSF determines the attractor 
dynamics to be synchronizable, while the actual simulation 
of a network of 6 oscillators never attains attractor dynamics. 
Nevertheless, this is an example of the importance of checking 
the practical applications of the MSF predictions, specially for 
high values of coupling strengths. 

When coupling is only introduced through the y variable 
(Fig. |^(b)) the oscillators behave as class II systems for low 
values of 77 , as expected. Nevertheless, for moderate values of 
the noise strength (3 < 77 < 5) the system shifts to a behav¬ 
ior similar to that of class III systems, synchronizing only for 
intermediate values of the coupling strength. Finally, when 
the noise is further increased (77 > 5, not shown here), the 
network is not able to reach the synchronized state for any 
value of the coupling, behaving as a class I system. This way, 
despite being a class II system (when coupled through the y 
variable), the addition of noise can make system behave dif¬ 
ferently. 

Finally, it is also worth analyzing how the combination of 
both layers increases the stability of the synchronized mani¬ 
fold. Figure [^(c) shows the synchronization error for dr = 
0.5, i.e. when the coupling is equally distributed among both 
X and Y layers. We can observe that, as suggested by the the¬ 
oretical predictions shown in Fig. the combined coupling 
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of the X and y variables enhances the network stability, so that 
the synchronization of the system is maintained even for high 
values of the noise strength. Obviously, only for small values 
of the coupling strength, here measured as u, unsynchronized 
motion is observed, as it is expected in class II systems. 

Note that the boundaries of the region where synchroniza¬ 
tion becomes robust show an excellent agreement with the 
zero crossings of the MSF. On the other hand, in the areas 
where synchronization is most robust against the presence of 
noise, the optimal v from that point of view is not always close 
to the minimum of the MSF. This is not surprising, consider¬ 
ing that we are using noise strengths of the order 77 = 5, which 
are certainly beyond any reasonable definition of infinitesi¬ 
mal perturbations, a necessary requirement for obtaining the 
MSF. The addition of finite perturbations, being a largely un¬ 
explored issue, are beyond of the scope of the MSF frame¬ 
work, and require specific numerical studies adapted to each 
particular topology, even if some general behaviour has been 
identified ms). 


IV. EXPERIMENTAL ANALYSIS 

To test the robustness of the theoretical predictions given 
by the MSF and the numerical simulations, we have imple¬ 
mented the electronic version of the Rossler system described 
in Eq. A schematic representation of the experimental de¬ 
sign is shown in Fig. It consists of an electronic array 
(FA), a personal computer (PC) and a multifunction data card 
(DAQ) composed by 12 analog to digital converters (ADCs) 
and 6 digital to analog converter (DACs). The ADCs are used 
for sampling the variable x of the oscillators, while the DACs 
generate six noisy signals that perturb the dynamics of each 
node separately. The EA comprises 6 Rossler-like electronic 
circuits forming a spiderweb network (blue nodes), with one 
central node and five peripheral nodes. Each node has two 
individual electronic couplers, one for the x variable and the 
second for the y variable, both controlled by a two digital po¬ 
tentiometers (XDCP), which are adjusted by a digital signal 
coming from ports PO.O - 1 (DO). PO.O is used to increase 
or decrease the resistance of the voltage divisor (a), and PO.l 
sets the value of the resistance (the resolution allowing for 100 
discretized steps). 

The noisy signals are designed in Labview, using the li¬ 
brary Gaussian White Noise VI 1^ that generates six differ¬ 
ent Gaussian-distributed pseudorandom sequences bounded 
between [—1,1]. These signals are digitally filtered by a third- 
order lowpass Butterworth filter 1^ with a cutoff frequency 
of 10 kHz. All the experimental process is controlled from a 
virtual interface developed in Labview 8.5. 

The experimental procedure is the following: first, a and 
T] are set to zero and then we introduce the six noisy signals 
and apply the factor gain ( 77 ). After a waiting time of 500 ms 
(roughly corresponding to P = 600 cycles of the autonomous 
systems), the signals corresponding to the variables of the 6 
circuits are acquired by the analog ports (AI 0; AI 1; ... ; AI 
5) and the synchronization error is calculated and stored in the 
PC. Noisy signals are injected by the digital converters (AO 0; 
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FIG. 5: (Color online) Experimental results. Synchronization errors 
as a function of the noise strength 77 and the coupling (normalized 
to a = dr 72 ). As in Fig. panels correspond to: (a) dr = 0, (b) 
dr = 1 and (c) dr = 0.5. Region colored in red indicates where 
the combination of noise and coupling strength leads to a loss of 
the complete synchronization. The network combining the coupling 
through the two variables x and 7 / (dr = 0.5) is the one showing 
better performance against noise perturbations. Left plots correspond 
to the MSF obtained for each value of dr. 


AO 1; ... ; AO 5) and this part of the process is repeated 100 
times (until the maximum value of a is reached). Finally, 77 is 
increased to the next value and a is swept again. The whole 
process is repeated 100 times until the maximum value of 77 is 
reached. 

Figure 1^ shows the experimental results for a configuration 
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identical to that of the numerical simulations shown in Fig. 

We observe that the qualitative agreement between nu¬ 
merics and experiment is excellent, in spite of unavoidable 
parameter mismatches in the experimental realization due to 
the tolerance of the electronic components (between 5% and 
10%). The parameter mismatch, together with the experi¬ 
mental noise, make the oscillators in the network not only 
slightly different from their mathematical definition, but also 
non-identical to one another. This way, we confirm experi¬ 
mentally the feasibility of using the MSF for evaluating how 
the coupling through multiple variables enhances the stability 
of the synchronous state of a network under realistic condi¬ 
tions. 


V. CONCLUSIONS 

We have seen how an adequate distribution of the coupling 
strength between the variables of a dynamical system leads 
to an enhancement of the stability of the synchronized mani¬ 
fold. In particular, we have shown that it is possible obtain a 
MSF that depends on the parameter dr accounting for the dis¬ 
tribution of strength, while maintaining the global coupling 
constant. Interestingly, we report the existence of an optimal 
value of dr indicating what is the most adequate amount of 
coupling to be considered at each coupling variable. The opti¬ 
mal value of dr is independent of the topology of the network, 
as long as we use the same coupling structure among all vari¬ 
ables. 

Using electronic circuits, we have also checked the robust¬ 
ness of the results when noise and parameter mismatch are 
considered, which confirms the theoretical predictions given 
by the parametric MSF and, in its turn, reveal that the require¬ 
ment of the oscillators to be identical can be relaxed. 

The proposed framework of decomposing the different di¬ 
mensions of the system (variables) in interconnected layers 
paves the way to the use of the multilayer networks tools 
12^ [23 to further analyze synchronization phenomena in 
multivariable coupled systems. Indeed, the current theoretical 
efforts in network theory to define and study complex struc¬ 
tures resulting from the interaction of networks, e.g. interde¬ 
pendent networks and multiplex networks among others, have 
made great progress in recent times, in showing new emergent 


phenomena with no counterpart in single (monolayer) com¬ 
plex networks ifTTI 12414^ . New developments in that direc¬ 
tion could be further extend the results we here present, either 
using, e.g., the insight developed within the hypemetwork for¬ 
malism uniini, or using new approaches for analyzing mul¬ 
tilayer networks GU. 

Indeed, our methodology has some limitations that must be 
further explored in the future. First of all, the fact that the 
coupling must have the same topological structure at all vari¬ 
ables is a strict constraint, since real systems may have differ¬ 
ent configurations depending on the coupling variable. More 
general, fully multilayer topologies could be considered by re¬ 
sorting to the hypernetwork formalism introduced in (TOl . and 
for greater generality one can use the method in Ref. C2. 
With this methodology, the case of A could be ad- 
dressed at the cost of introducing some more complexity to 
the problem. Nevertheless, it would be of great interest, since 
it would raise new questions such as what the adequate com¬ 
bination of topologies would be given a specific distribution 
of weights dr. Second, since the parametric MSF depends on 
the dynamical system implemented in the network, we can not 
guarantee the existence of an optimal balance of the distribu¬ 
tion of coupling between layers in other dynamical systems, 
at least until their corresponding MSF have been analyzed. 

VI. ACKNOWLEDGEMENTS 

Authors acknowledge D. Papo and P.L. del Barrio for fruit¬ 
ful conversations. Support from MINECO through projects 
FIS2011-25167, FIS2012-38266 and FIS2013-41057-P is 
also acknowledged. AA and JGG acknowledge support from 
the EC EET-Proactive Project PLEXMATH (grant 317614) 
and MULTIPLEX (grant 317532). JGG acknowledges sup¬ 
port from MINECO through the Ramon y Cajal program, 
the Comunidad de Aragon (Grupo EENOL) and the Brazil¬ 
ian CNPq through the PVE project of the Ciencia Sem 
Eronteiras program. AA acknowledges ICREA Academia 
and the James S. McDonnell Eoundation. R.S.E. ac¬ 
knowledges Universidad de Guadalajara, CULagos (Mexico) 
for financial support (OP/PIEI-2013-14MSU0010Z-17-04, 
PROINPEP-RG/005/2014, UDG-CONACyT/IOlO/163/2014) 
and CONACyT (Becas Mixtas MZO2015/290842). 


[1] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez and D.-U. 
Hwang, Phys. Rep. 424,175-308 (2006). 

[2] A. Arenas, A. Diaz-Guilera, J. Kurths, Y. Moreno and C. Zhou, 
Phys. Rep. 469, 93-153 (2008). 

[3] A. Barrat, M. Barthelemy and A. Vespignani. Dynamical Pro¬ 
cesses on Complex Networks, Cambridge University Press 
(2008). 

[4] L.M. Pecora and T.L. Carroll, Phys. Rev. Lett. 80, 2109 (1998). 

[5] A.E. Motter, C. Zhou and J. Kurths, Phys. Rev. E 71, 016116 
(2005). 

[ 6 ] M. Chavez, D.-U. Hwang, A. Amann, H. G. Hentschel and S. 
Boccaletti, Phys. Rev. Lett. 94, 218701 (2005). 


[7] A. E. Motter, C.S. Zhou and J. Kurths, Europhys. Lett. 69, 334- 
340 (2005). 

[ 8 ] T. Nishikawa and A.E. Motter, Phys. Rev. E 73, 065106 (2006). 

[9] L. Donetti, P. I. Hurtado and M. A. Munoz, Phys. Rev. Lett. 95, 
188701 (2005). 

[10] F. Sorrentino, New J. Phys. 14, 033035 (2012). 

[ 11 ] A network with a connectivity matrix A = {cHj} (cHj = 1 
if nodes i and j are connected, and 0 otherwise) can be rep¬ 
resented by the Laplacian matrix C = {hj}, where kj = 

~ denoting the Kronecker delta). 

[12] D. Irving and F. Sorrentino, Phys. Rev. E 86, 056102 (2012). 

[13] A.N. Pisarchik, R. Jaimes-Reategui, J.R. Villalobos-Salazar, 



8 


J.H. Garcia-Lopez and S. Boccaletti, Phys. Rev. Lett. 96, 
244102 (2006). 

[14] J. Bragard, S. Boccaletti and H. Mancini, Phys. Rev. Lett. 91, 
064103 (2003). 

[15] D.-U. Hwang, M. Chavez, A. Amann and S. Boccaletti, Phys. 
Rev. Lett. 94, 138701 (2005). 

[16] L. Huang, Q. Chen, Y.-C. Lai and L.M. Pecora, Phys. Rev. E 
80, 036204 (2009). 

[17] J. Aguirre, R. Sevilla-Escoboza, R. Gutierrez, D. Papo and J. 
M. Buldii, Phys. Rev. Lett. 112, 248701 (2014). 

[18] R. Sevilla-Escoboza, J. M. Buldii, A. N. Pisarchik, S. Boc¬ 
caletti, and R. Gutierrez, Phys. Rev. E 91, 032902 (2015). 

[19] R. Gutierrez, E. del-Pozo and S. Boccaletti, PLoS ONE 6, 
e20236 (2011). 

[20] National Instrument, Gaussian White Noise VI 

http://zone.ni.com/reference/en-XX/help/ 
3713 61J-01/Ivanls/gaussian_white_noise/ 


[21] National Instrument, Filter http://zone.ni. 
com/reference/en-XX/help/3713 61J-01/ 
lvexpress/signal_filter/ 

[22] M. Kivela, A. Arenas, M. Barthelemy, J.P. Gleeson, Y. Moreno, 
M.A. Porter, Multilayer Networks J. Complex Netw. 2, 203 
(2014) 

[23] S. Boccaletti, G. Bianconi, R. Criado, C.I. del Genio, J. Gomez- 
Gardenes, M. Romance, I. Sendina-Nadal, Z. Wangk, M. Zanin, 
Phys. Rep. 544, 1 (2014). 

[24] J. Gao, S.V. Buldyrev, H.E. Stanley and S. Havlin, Nat. Phys. 8 , 
40 (2012). 

[25] J. Aguirre, D. Papo and J. M. Buldii, Nat. Phys. 9, 230 (2013). 

[26] E. Radichi and A. Arenas, Nat. Phys., 9 717 (2013). 

[27] M. De Domenico, A. Sole-ribalta, E. Cozzo, M. Kivela, Y. 
Moreno, M.A. Porter, S. Gomez and A. Arenas, Phys. Rev. X 
3, 041022 (2013). 


